clear; close all; clc;

% 模型参数
alpha = 4.5;       % (α, σ, μ) = (4.5, 0.14, 0.001)
sigma = 0.14;
mu    = 0.001;

% 耦合方程中 phi 的更新系数
epsilon = 1;       % ε = 1

N_total = 50000;   % 总迭代次数
N_cut   = 10000;   % 去除瞬态的步数

%% (a) 扫描 k，固定 φ₀ = -0.3
phi0_fixed = -0.3;
k_range = linspace(-0.4, 0.4, 101);  
E_a = zeros(1, length(k_range));

for i = 1:length(k_range)
    k_val = k_range(i);
    [~, ~, E_a(i)] = Rulkov(alpha, sigma, mu, k_val, phi0_fixed, N_total, epsilon, N_cut);
end

%% (b) 扫描 φ₀，固定 k = 0.1
k_fixed = 0.1;
phi0_range = linspace(-2, 2, 101);  
E_b = zeros(1, length(phi0_range));

for i = 1:length(phi0_range)
    phi_val = phi0_range(i);
    [~, ~, E_b(i)] = Rulkov(alpha, sigma, mu, k_fixed, phi_val, N_total, epsilon, N_cut);
end

%% 绘图
figure;

% (a) 子图：E 随 k 变化，固定 φ₀ = -0.3
subplot(1,2,1);
plot(k_range, E_a, 'LineWidth',1 , 'Color',[0.7,0.3,0.0]);
xlabel('k');
ylabel('E');
title('(a) \phi_0 = -0.3');

% (b) 子图：E 随 φ₀ 变化，固定 k = 0.1
subplot(1,2,2);
plot(phi0_range, E_b, 'LineWidth',1 , 'Color',[0.0,0.6,0.0]);
xlabel('\phi_0');
ylabel('E');
title('(b) k = 0.1');
